Variational study of hard-core bosons in a 2— D optical lattice using Projected 

Entangled Pair States (PEPS) 
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We have studied the system of hard-core bosons on a 2-D optical lattice using a variational 
algorithm based on projected entangled-pair states (PEPS). We have investigated the ground state 
properties of the system as well as the responses of the system to sudden changes in the parameters. 
We have compared our results to mean field results based on a Gutz wilier ansatz. 

PACS numbers: 03.75.Lm, 02.70.-c, 75.40.Mg, 03.67.-a 



I. INTRODUCTION 

Systems of interacting bosons in optical lattices have 
attracted a lot of interest in the last few years due to 
recent experimental achievements [l|, H, H, H, [E @1 • In 
these systems, atoms are trapped by the combination of 
a periodic and a harmonic potential created by coun- 
terpropagating laser-beams. These systems of trapped 
atoms resemble a crystal in the sense that atoms are lo- 
calized at periodic locations. The theoretical study of the 
static and dynamic properties of these systems is quite 
challenging. An exact solution can only be obtained in 
one dimension in the so-called Tonks-Girardeau limit via 
fermionization 0, 0] ■ Outside this limit and in higher di- 
mensions, approximate methods have to be used. The 
most common of those, the mean field approximation 
based on a Gutzwiller ansatz Q , is unfortunately known 
to give imprecise predictions for correlations. The most 
powerful numerical methods such as the Density Ma- 
trix Renormalization Group (DMRG) @ and Quantum 
Monte Carlo (QMC) are also of limited use: DMRG is 
mainly restricted to 1-D systems and QMC suffers from 
the sign-problem as time-evolutions are investigated. 

In this paper, we apply the algorithm introduced 
in [lol ]. This algorithm is a variational method within 
a class of states termed Projected Entangled-Pair States 
(PEPS). It has been proven to work well for the Heisen- 
berg antiferromagnet flfjj and the frustrated Shastry- 
Sutherland model The system we focus on is the 

system of hard-core bosons in a 2-D optical lattice. This 
model captures the essential physics of bosons in opti- 
cal lattices [l2| . We use the algorithm to determine the 
ground state properties of the system and to study the 
responses of the system to sudden changes in the pa- 
rameters. We compare our results to mean field results 
based on the Gutzwiller ansatz. We find that the PEPS 
and the Gutzwiller ansatz deviate clearly in the predic- 
tion of the ground state momentum distribution and the 
time-evolution of the condensed fraction of the particles. 

The paper is outlined as follows: In section [II] we give 
a brief overview of the algorithm. We specify the model 
we want to investigate in section Hill and explain the way 
the algorithm is applied to this model in section HVl The 



results of the numerical calculations are presented in sec- 
tions [V] and IVIl We conclude with the discussion of the 
the performance and the stability of the algorithm in sec- 
tion ELD 



II. THE ALGORITHM 

Although the algorithm has already been outlined 
in [lfj , we reiterate it here in detail - specialized for our 
purpose. The algorithm is a variational method with re- 
spect to the class of PEPS. These states have been found 
to be adequate for representing the ground state of nu- 
merous many-body systems. Also, these states are favor- 
able to variational calculations because they possess an 
internal refinement parameter, the virtual dimension D, 
that allows to control the precision of the calculation. 
While D = 1 specializes the PEPS to a product state, 
the choice D = d M (with M being the total number of 
lattice sites and d the dimension of one subsystem) en- 
larges the space of PEPS to the complete Hilbert-space 
of the system. The purpose of the algorithm is - in our 
case - to simulate the time-evolution of a system within 
the subset of PEPS with a fixed D. This means that 
after each time-evolution step the state of the system 
is approximated by the "nearest" PEPS with virtual di- 
mension D. The key element of the algorithm is thus the 
optimization of the parameters of a PEPS such that its 
distance to a given state is minimized. 

The manner in which the optimization is performed 
is closely related to the structure of PEPS. A PEPS is 
a state with coefficients that are contractions of tensors 
according to a certain scheme. Thereby, each tensor is 
associated with a physical subsystem. The contraction- 
scheme mimics the underlying lattice structure. Each 
tensor possesses one physical index with dimension equal 
to the physical dimension d of a subsystem and a certain 
number of virtual indices with dimension D. The number 
of virtual indices is equal to the number of bonds that em- 
anate from the lattice-site the tensor is associated with. 
For example, in a rectangular lattice the number of vir- 
tual indices is 4 (except at the borders) - related to the 
left, right, upper and lower bond respectively. The tensor 
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FIG. 1: Structure of the coefficient related to the state 
I fen, .., fei4 ) in the PEPS I^a). The bonds represent the 
indices of the tensors [Ai] k that are contracted. 



FIG. 2: Structure of the contractions in { \E" a \ )■ In this 
scheme, the first and last rows can be interpreted as MPS 
I Ui } and ( U4 I and the rows in between as MPO U2 and U3. 
The contraction of all tensors is then equal to ( f/4 If/3 (72 1 Ui ). 



associated with site i is 

M Irud 

with physical index k and virtual indices I, r, u and d. 
The coefficients of the PEPS are then formed by joining 
the tensors in such a way that all indices related to same 
bonds are contracted. This is illustrated in fig. Q] for 
the special case of a 4 x 4 square lattice. Assuming this 
contraction of tensors is performed by the function •?"(•), 
the resulting PEPS can be written as 

d 

fel,...,feM=l 

The aim of the algorithm is to optimize the tensors Ai 
such that the distance between the PEPS | ) and a 
given state tends to a minimum. We assume that the 
given state is a PEPS \^b) with virtual dimension Db 
and tensors Bi. This is no loss of generality since ev- 
ery state can be written as a PEPS. The function to be 
minimized is then 

K(A U ...,A M ) = |||* A )-|* B )|| 2 . 

This function is non-convex with respect to all parame- 
ters {Ai, A M }. However, due to the special structure 
of PEPS, it is quadratic in the parameters Ai associated 
with one lattice-site i. Because of this, the optimal pa- 
rameters Ai can simply be found by solving a system of 
linear equations. The concept of the algorithm is to do 
this one-site optimization site-by-site until convergence 
is reached. 

The challenge that remains is to calculate the co- 
efficient matrix and the inhomogenity of the linear 
equations-system. In principle, this is done by contract- 
ing all indices in the expressions for the scalar-products 



( Of a 1 <3>a ) and ( ^a \^b) except those connecting to Aj. 
By interpreting the tensor A4 as a dD 4 -dimensional vec- 
tor Ai, these scalar-products can be written as 

= AtMAi (1) 
(*a\*b) = A\Wi. (2) 

Since 

K= (*B|*B) + <*A|*A)-2Jfe(* A |* B ), 

the minimum is attained as 

A/jij = Wi. 

The obstacle, however, is that the numerical calculation 
of the coefficient matrix Mi and the inhomogenity VVj re- 
quires a number of operations that scales exponentially 
with the number of subsystems M. This will make the 
algorithm non-efficient as the system grows larger. Be- 
cause of this, an approximate method has to be used to 
calculate Mi and VVj . 

The approximate method suggested in 10] is based 
on matrix product states (MPS) and matrix product op- 
erators (MPO). To see how MPS and MPO implicitly 
appear in the problem of calculating Mi and VVj, we take 
a closer look at the structure of the contractions in the 
scalar-products ( 4" a \ ^ A ) and ( '5 a | $b )■ Thereby, we 
focus on a L x L square lattice in the following. We start 
with the study of ( *S?a | ^a)- For this, we single out a 
specific site j and define the D 2 x D 2 x D 2 x D 2 -tensor 

d 

r -I {uu'){dd') _ \ - [ A*l k { A ~\ k 
L J-l(M')Or') — L jllrudl Hl'r'u'd'' 
k=l 
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In this definition, the symbols {IV), (rr'), (uu') and (dd!) 
indicate composite indices. We may interpret the 4 in- 
dices of this tensor as being related to the 4 bonds em- 
anating from site j in the lattice. Then, ( ^a | ^ A ) is 
formed by joining all tensors Ej in such a way that all 
indices related to same bonds are contracted - as in the 
case of the coefficients of PEPS. These contractions have 
a rectangular structure, as depicted in fig. [21 In terms of 
the function .F(-), the scalar-product reads 

(*a|*a) =J 7 (E 1 ,...,E m ). 

The main idea of the approximate algorithm is to in- 
terpret the first and the last row in this contraction- 
structure as MPS and the rows in between as MPO. The 
horizontal indices thereby form the virtual indices and 
the vertical indices are the physical indices. Thus, the 
MPS and MPO have both virtual dimension and physi- 
cal dimension equal to D 2 . Explicitly written, the MPS 
read 

| CM = E tr([E 11 f 1 ---[E 1L f L )\d 1 ,...J L ) 

d 1 ,...,d L = l 

(U L \ = tr([E L1 ]* ll ...[E LI ]**yu 1 ,...,u L \ 
and the MPO at row r is 

U r = tv([E rl ^-.-[E rL ]^ L )x 

ui,...,ii=l 
dt,...,d L = l 

x \u 1 ,...,u L ){d ll ...,d L \. 

In terms of these MPS and MPO, the scalar-product is 
a product of MPO and MPS: 

= {Ul\U L -i---U 2 \Ux) 

The evaluation of this expression is, of course, in- 
tractable. With each multiplication of a MPO with a 
MPS, the virtual dimension increases by a factor of D 2 . 
Thus, after L multiplications, the virtual dimension 
is D 2L - which is exponential in the number of rows. 
The expression, however, reminds of the time-evolution 
of a MPS [H, [H Gj|. There, each multiplication with 
a MPO corresponds to one evolution-step. The prob- 
lem of the exponential increase of the virtual dimension 
is circumvented by restricting the evolution to the sub- 
space of MPS with a certain virtual dimension D. This 
means that after each evolution-step the resulting MPS 
is approximated by the "nearest" MPS with virtual di- 
mension D. This approximation can be done efficiently, 
as shown in [l3j]. In this way, also a\^ a) can be 
calculated efficiently: first, the MPS | U 2 ) is formed by 
multiplying the MPS | t/i ) with MPO U 2 . The MPS | U 2 ) 
is then approximated by | U 2 ) with virtual dimension D. 



In this fashion the procedure is continued until | Ul-i) is 
obtained. The scalar-product ( <P a \ ^ A ) is then simply 

= (U l \Ul-i). 

The calculation of the coefficient matrix Mi is closely re- 
lated to the calculation of ( ^ a \ ^ A ) ■ Ml relies on the 
contraction of all but one of the tensors Ej according to 
the same scheme as before. The one tensor that has to 
be omitted is E\ - the tensor related to site i. Assuming 
this contraction is performed by the function Qi{-), Mi 
can be written as 

Wilt r , l:' U ' d ' =Qi(Ei,...,E M ) l '/' u ' d ' 8%. 

L '■llrud k' c v 11 ' lvl J Irud K 

If we join the indices (klrud) and (k'l'r'u'd 1 ), we obtain 
the dD 4 x dD 4 -matrix that fulfills equation (TT|) . To evalu- 
ate Gi(-) efficiently, we proceed in the same way as before 
by interpreting the rows in the contraction-structure as 
MPS and MPO. First, we join all rows that lie above 
site i by multiplying the topmost MPS \ U\) with subja- 
cent MPO and reducing the dimension after each multi- 
plication to D. Then, we join all rows lying below i by 
multiplying ( Ul | with adjacent MPO and reducing the 
dimension as well. We end up with two MPS of virtual 
dimension D - which we can contract efficiently with all 
but one of the tensors Ej lying in the row of site i. 

The scalar-product ( a \^b) and the inhomogenity 
Wi are calculated in an efficient way following the same 
ideas. First, the DDb x DDb X DDb x CDs-tensors 

d 

r i(uu')(dd') _ ^ \ A*] k \T3] k 

rJ'J (H')(rr') V 3\lrudV D nVr>u'd'- 

fc=l 

are defined. The scalar-product (^a \^b) is then ob- 
tained by contracting all tensors Fj according to the pre- 
vious scheme - which is performed by the function T{-): 

(V A \V B )=F(F 1 ,...,F M ) 

The inhomogenity Wi relies on the contraction of all but 
one of the tensors Fj, namely the function Q%y), in the 
sense that 

D 

i W i\lrud= ll &( F li-> F M)J u d l B i\l'r'u>d>- 
V v' u' d' — 1 

Joining all indices (klrud) in the resulting tensor leads 
to the vector of length dD 4 that fulfills equation (|2|). 
Thus, both the scalar-product ( ^a \ ) and the in- 
homogenity Wi are directly related to the expressions 
!F{Fi, F M ) and Gi{F\, Fm)- These expressions, 
however, can be evaluated efficiently using the approx- 
imate method from before. 

Summing up, we have an algorithm that allows the 
efficient reduction of the virtual dimension of a PEPS - 
and thus the efficient simulation of a time-evolution step 
within the subset of PEPS. 
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III. THE MODEL: HARD-CORE BOSONS IN A 
2 D OPTICAL LATTICE 

This algorithm we use to study a system of bosons 
in a 2-D optical lattice of size L x L. This system is 
characterized by the Bose-Hubbard Hamiltonian 

H = - J E { a l a j + h.c) + y hjlfij - 1) + E Vjfk, 

<i,j> i i 

where a\ and are the creation and annihilation opera- 
tors on site i and hi — ajeij is the number operator. This 
Hamiltonian describes the interplay between the kinetic 
energy due to the next-neighbor hopping with ampli- 
tude J and the repulsive on-site interaction U of the par- 
ticles. The last term in the Hamiltonian models the har- 
monic confinement of magnitude Vi = Vo(i — «o) 2 - Since 
the total number of particles N — hi is a symmetry of 
the Hamiltonian, the ground-state will have a fixed num- 
ber of particles N. We choose this number by appending 
the term —fiN to the Hamiltonian and tuning the chem- 
ical potential /i. The variation of the ratio U/J drives a 
phase-transition between the Mott-insulating and the su- 
perfluid phase, characterized by localized and delocalized 
particles respectively [l6j . Experimentally, the variation 
of U/J can be realized by tuning the depth of the optical 
lattice [13, [H| . The quantity that is typically measured 
is the momentum distribution. The is done by letting 
the atomic gas expand and measuring the density distri- 
bution of the expanded cloud. Thus, we will be mainly 
interested here in the (quasi)-momentum distribution 

r.s 

of the particles. 

In the following, we focus on the limit of a hard-core 
interaction, U/J — > oo. In this limit, two particles are 
prevented from occupying a single site. This limit is es- 
pecially interesting in one dimension where the particles 
form the so-called Tonks-Girardeau gas 0, 0] • The par- 
ticles in this gas are strongly correlated - which leads to 
algebraically decaying correlation functions. In two di- 
mensions, the model was studied in detail in In the 
hard-core limit, the Bose-Hubbard model is equivalent 
to a spin-system with XX-interactions described by the 
Hamiltonian 

h = ~ E + + \ E (v* - 

<i,j> i 

Here, Ox \ ay and denote the Pauli-operators act- 
ing on site i. This Hamiltonian has the structure we 
can simulate with the algorithm: it describes L 2 physical 
systems of dimension d — 2 on a L x L-square lattice. 



IV. APPLICATION OF THE ALGORITHM TO 
HARD-CORE BOSONS IN A 2-D LATTICE 

The principle of simulating a time-evolution step ac- 
cording to XX-Hamiltonian is as follows: first, a PEPS 
| *f? A ) with physical dimension d = 2 and virtual dimen- 
sion D is chosen as a starting state. This state is evolved 
by the time-evolution operator U = e~ lHSt (we assume 
H = 1) to yield another PEPS | VPs ) with increased vir- 
tual dimension Db- 

\*b) = U\*° a ) 

The virtual dimension of this state is then reduced to D 
by applying the algorithm of the previous section. This 
means, a new PEPS | ) with virtual dimension D 
is calculated that has minimal distance to | ^>b ) ■ This 
new PEPS is then the starting state for the next time- 
evolution step. 

The operator U, however, increases the virtual dimen- 
sion of a PEPS by a factor that scales exponentially 
with L. This is why it is more convenient to approxi- 
mate U by an operator that increases the virtual dimen- 
sion merely by a constant factor 77. This is done by means 
of a Trotter-approximation: first, the interaction-terms 
are classified in horizontal and vertical according to their 
orientation and in even and odd depending on whether 
the interaction is between even-odd or odd-even rows 
(or columns). The Hamiltonian can then be decomposed 
into a horizontal-even, a horizontal- odd, a vertical-even 
and a vertical-odd part: 

H = H] le + Hho + H ve + H vo 

The single-particle operators of the Hamiltonian can sim- 
ply be incorporated in one of the four parts. Using the 
Trotter-approximation, the time-evolution operator U 
can be written as a product of four evolution-operators: 

jj = e -iHSt e -iH ha St e -iH ho 8t e ~iHv C St e ~tH vo 8t /g\ 

Since each of the four parts of the Hamiltonian consists 
of a sum of commuting terms, each evolution-operator 
equals a product of two-particle operators 

acting on neighboring sites i and j. These two-particle 
operators have a Schmidt-decomposition consisting of 4 
terms: 

4 

w y = E u i ® V J 
P =i 

One such two-particle operator Wij applied to the 
PEPS I ^ ) modifies the tensors A® and A® associated 
with sites i and j as follows: assuming the sites i and j 
are horizontal neighbors, A® has to be replaced by 

k 2 k ' 

\. Bj \ l(rp)ud = E \- U i\ W Irud 
k'=l 
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FIG. 3: Energy as a function of time for the imaginary time- 
evolution of the system of hard-core bosons on a 4 x 4-lattice. 
The evolutions are performed sequentially with PEPS of vir- 
tual dimension D — 2, D — 3, D — 4 and D = 5. The times 
at which D is increased are indicated by vertical lines. For 
comparison, the exact ground state-energy, the exact imagi- 
nary time-evolution and the energy of the optimal Gutzwiller 
ansatz are included. 



and A® becomes 



k 

(lp)rud 



E 



"J \k 



[A? 



k! 

J J Irud' 



These new tensors have a joint index related to the bond 
between sites i and j. This joint index is composed of 
the original index of dimension D and the index p of 
dimension 4 that enumerates the terms in the Schmidt- 
decomposition. Thus, the effect of the two-particle oper- 
ator Wij is to increase the virtual dimension of the bond 
between sites i and j by a factor of 4. Consequently, 
e -iH he 8t anc j e -tH ho st mcrease t ne dimension of every 
second horizontal bond by a factor of 4; e ~ lH i^ St anc l 
e -%n vo U fag same for every second vertical bond. By 
applying all four evolution operators consecutively, we 
have found an approximate form of the time-evolution 
operator U that - when applied to a PEPS - yields an- 
other PEPS with a virtual dimension multiplied by a 
constant factor r\ = 4. 

Even though the principle of simulating a time- 
evolution step has been recited now, the implementation 
in this form is numerically expensive. This is why we ap- 
pend some notes about how to make the simulation more 
efficient: 

1.- Partitioning of the evolution: The number of required 
numerical operations decreases significantly as one time- 
evolution step is partitioned into 4 substeps: first the 
state | ) is evolved by e - tH voSt on iy anc i th e dimen- 
sion of the increased bonds is reduced back to D. Next, 



evolutions according to e - lH "^ 5t : e -iH ho 5t an( ^ e -iH he 8t 
follow. Even though the partitioning increases the num- 
ber of evolution-steps by a factor of 4, the number of 
multiplications in one evolution-step decreases by a fac- 
tor of rf = 64. 

2. - Optimization of the contraction order: Most critical 
for the efficiency of the numerical simulation is the order 
in which the contractions are performed. We have op- 
timized the order in such a way that the scaling of the 
number of multiplications with the virtual dimension D is 
minimal. For this, we assume that the dimension D that 
tunes the accuracy of the approximate calculation of M% 
and Wj is proportional to D 2 , i.e. D = kD 2 . The num- 
ber of required multiplications is then of order k 2 D 12 L 2 
and the required memory scales as dr/K 2 D 8 . 

3. - Optimization of the starting state: The number of 
sweeps required to reach convergence depends on the 
choice of the starting state for the optimization. The 
idea for finding a good starting state is to reduce the 
bonds with increased virtual dimension rjD by means of 
a Schmidt-decomposition. This is done as follows: as- 
suming the bond is between the horizontal neighboring 
sites i and j, the contraction of the tensors associated 
with these sites, Bi and Bj, along the bond i—j forms 
the tensor 



[M 



k k' 
l 3\lud r'u'd' 



Or) 

E 

P =i 



B, 



1 M fc ' 



l pud I 



By joining the indices (Mud) and {k'r'u'i 
can be interpreted as a AD 3 x eLD 3 -matrix. 
decomposition of this matrix is 



pr'u' d' ' 

I'), this tensor 
The Schmidt- 



dD J 



= y c,,a: • 



with the Schmidt-coefficients c p (c p > 0) and correspond- 
ing matrices Al and A? . We can relate these matrices to 
a new pair of tensors A® and A® associated with sites i 
and j: 



[Al 



,k 

I l pud 
k 

prud 



k 

lud 
k 

rud 



The virtual dimension of these new tensors related to 
the bond between sites i and j is equal to the num- 
ber of terms in the Schmidt-decomposition. Since these 
terms are weighted with the Schmidt-coefficients c p , it 
is justified to keep only the D terms with coefficients 
of largest magnitude. Then, the contraction of the ten- 
sors A® and A® along the bond i—j with dimension D 
yields a good approximation to the true value Mi 3 ■: 



[M 



k' 

1 3 J lud r'u'd' 



D 

E 

P =i 



[ A i] Ipud [ A °j] pr'u' d' ' 



This method applied to all bonds with increased dimen- 
sion provides us with the starting state for the optimiza- 
tion. 
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FIG. 4: Energy as a function of time for the imaginary time- 
evolution of the system of hard-core bosons on a 11 x 11- 
lattice. The evolutions are performed sequentially with PEPS 
of virtual dimension D = 2, D — 3, D = 4 and D = 5. The 
times at which D is increased are indicated by vertical lines. 
For comparison, the energy of the optimal Gutzwiller ansatz 
is included. 



V. GROUND STATE PROPERTIES 

In the following, we study the ground-state properties 
of the system of hard-core bosons for lattice-sizes 4x4 
and 11 x 11. We calculate the ground-state by means of 
an imaginary time-evolution which we can simulate with 
the method from before. 

We first focus on the 4 x 4-lattice for which we can cal- 
culate the ground-state exactly and are able to estimate 
the precision of the algorithm by comparison with exact 
results. In fig. [31 the energy is plotted as the system 
undergoes the imaginary time-evolution. We thereby as- 
sume a time-step St = — z0.03. We choose the magnitude 
of the harmonic confinement (in units of the tunneling- 
constant) Vo/J — 36. In addition, we tune the chem- 
ical potential to fi/J = 3.4 such that the ground state 
has particle-number N = 4. With this configuration, 
we perform the imaginary time-evolution both exactly 
and variationally with PEPS. As a starting state we 
take a product state that represents a Mott-like distri- 
bution with 4 particles arranged in the center of the trap 
and none elsewhere. The variational calculation is per- 
formed with D = 2 first until convergence is reached; 
then, evolutions with D — 3, D — 4 and D = 5 fol- 
low. At the end, a state is obtained that is very close 
to the state obtained by exact evolution. The difference 
in energy is \E D=5 — E exact \ ^ 6.4614 • 10~ 5 J. For com- 
parison, also the exact ground-state energy obtained by 
an eigenvalue-calculation and the energy of the optimal 
Gutzwiller ansatz are included in fig. [3J The difference 



FIG. 5: (Quasi)-momentum distribution of the particles in 
the ground state of a 11 x 11-lattice. Plotted are results of 
the variational calculations with PEPS of dimension D = 5 
and with the Gutzwiller ansatz. From the inset, the density 
of the particles can be gathered. 

between the exact result and the results of the imaginary 
time-evolution is due to the Trotter-error and is of order 
0(St 2 ). The energy of the optimal Gutzwiller- Ansatz is 
well seperated from the exact ground-state energy and 
the results of the imaginary time-evolution. 

In fig. [U the energy as a function of time is plotted 
for the imaginary time-evolution on the 11 x 11-lattice. 
Again, a time-step 8t = — i0.03 is assumed for the evolu- 
tion. The other parameters are set as follows: the ratio 
between harmonic confinement and the tunneling con- 
stant is chosen as Vq/J = 100 and the chemical potential 
is tuned to [ij J — 3.8 such that the total number of parti- 
cles N is 14. The starting state for the imaginary time- 
evolution is, similar to before, a Mott-like distribution 
with 14 particles arranged in the center of the trap. This 
state is evolved within the subset of PEPS with D = 2, 
D = 3, D — 4 and D — 5. As can be gathered from 
the plot, this evolution shows a definite convergence. In 
addition, the energy of the final PEPS lies well below the 
energy of the optimal Gutzwiller ansatz. 

The difference between the PEPS and the Gutzwiller 
ansatz becomes more evident as one studies the momen- 
tum distribution of the particles. The diagonal slice of 
the (quasi)-momentum distribution is shown in fig. EI As 
can be seen, there is a clear difference between the mo- 
mentum distribution derived from the PEPS and the one 
from the Gutzwiller ansatz. In contrast, the PEPS and 
the Gutzwiller ansatz produce a very similar density pro- 
file (see inset). The acceptability of the Gutzwiller ansatz 
is due to the inhomegenity of the system: the different 
average particle number at each site is the cause for the 
correlations between different sites. These correlations 
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FIG. 6: Time evolution of the condensate density after a 
sudden change of the magnitude of the trapping potential 
from Vo/J = 100 to Vo/J = 64. As a starting state, we 
use the Gutzwiller-approximation of the ground state. The 
evolution is performed on the basis of the Gutzwiller ansatz 
and PEPS with D = 2, D = 3 and D = 4. From the inset, 
the overlap between the PEPS with D — 2 and D = 3 (solid 
line) and the PEPS with D = 3 and D — 4 (dashed line) can 
be gathered. 



are, in many cases, good approximations. In contrast, 
the average particle number is constant in homogeneous 
systems - which leads to correlations that are constant. 
Thus, the Gutzwiller ansatz is expected to be less ap- 
propriate for the study of correlations of homogeneous 
systems. 



VI. DYNAMICS OF THE SYSTEM 

We now focus on the study of dynamic properties of 
hard-core bosons on a lattice of size 11 x 11. We in- 
vestigate the responses of this system to sudden changes 
in the parameters and compare our numerical results to 
the results obtained by the Gutzwiller ansatz. The prop- 
erty we are interested in is the fraction of particles that 
are condensed. For interacting and finite systems, this 
property is measured best by the condensate density p 
which is defined as largest eigenvalue of the correlation- 
matrix (a|aj). 

First, we study the time evolution of the condensate 
density after a sudden change of the trapping poten- 
tial. We start with a Gutzwiller-approximation of the 
ground state in case of a trapping potential of magni- 
tude Vq/J — 100. The chemical potential we tune to 
fj,/J = 3.8 to achieve an average particle-number of 
(N) = 14. This state we expose to a trapping potential 
of magnitude Vq/J = 64 and calculate the evolution of 



FIG. 7: Time evolution of the condensate density after a 
sudden shift of the center of the trap by one site in x— and 
y-direction. Starting state is the Gutzwiller-approximation 
of the ground state. The evolution is performed using the 
Gutzwiller ansatz and PEPS with D = 2, D = 3 and D = 4. 
The inset shows the overlap between the PEPS with D = 2 
and D = 3 (solid line) and D = 3 and D — 4 (dashed line). 



the condensate density using the Gutzwiller ansatz and 
PEPS with D = 2, D = 3 and D = 4. We thereby 
assume a time-step St = 0.03. To assure that our re- 
sults are accurate, we proceed as follows: first, we per- 
form the simulation using PEPS with D = 2 and D = 3 
until the overlap between these two states falls below a 
certain value. Then, we continue the simulation using 
PEPS with D = 3 and D = 4 as long as the overlap 
between these two states is close to 1. The results of this 
calculation can be gathered from fig. [6] What can be ob- 
served is that the results obtained from using PEPS are 
qualitatively very different from the result based on the 
Gutzwiller ansatz. The inset in fig. [6] shows the overlap 
of the D — 2 with the D = 3-PEPS and the D — 3 with 
the D = 4-PEPS. 

In fig. [7J the time-evolution of the condensate density 
after a sudden shift of the trapping potential is plotted. 
As a starting state, again the Gutzwiller-approximation 
of the ground state in a trap of magnitude Vo/J = 100 
is used. This state is evolved with respect to a trap- 
ping potential that is shifted by one lattice-site in x— 
and y-direction. We assume a time-step St — 0.03 and 
tune the chemical potential to fi/J= 3.8. As before, we 
perform the simulation successively with D — 2, D = 3 
and D = 4 and judge the accuracy of the results by 
monitoring the overlap between PEPS with different Ds. 
From the plot, it can be gathered that the evolution of 
the condensate density based on the Gutzwiller ansatz is 
qualitatively again very different from the evolution ob- 
tained from using PEPS. The evolution obtained from 
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FIG. 8: Time evolution of the condensate density starting 
from a Mott-distribution with 14-particles arranged in the 
center of the trap. The magnitude of the trapping potential 
is Vo/J= 100. For the evolution, the Gutzwiller ansatz and 
PEPS with D = 2, D = 3 and D = 4 are used. The inset 
shows the overlap between the D — 2 and D = 3-PEPS (solid 
line) and the D = 3 and D = 4-PEPS (dashed line). 



FIG. 9: Distance K between the time-evolved state and the 
state with reduced virtual dimension. The virtual dimensions 
D — 2, D = 3 and D — 4 are included. The distance is 
plotted for the evolution of a Mott-distribution with N = 14, 
as explained in fig. [8] From the inset, the deviation of the 
particle number from the value 14 can be gathered. 



using PEPS shows a definite damping. The shift of the 
trap thus provokes a destruction of the condensate. The 
evolution based on the Gutzwiller ansatz doesn't show 
this feature. 

As a contrary example, we study the evolution of a 
Mott-distribution with 14 particles arranged in the center 
of the trap. We assume Vo/J = 100, fj,/J = 3.8 and 
6t = 0.03. We perform the simulation in the same way 
as before with D = 2, D = 3 and D = 4. In fig. H the 
time evolution of the condensate density is plotted. It 
can be observed that there is a definite increase in the 
condensate fraction. The Gutzwiller ansatz is in contrast 
to this result since it predicts that the condensate density 
remains constant. 



VII. ACCURACY AND PERFORMANCE OF 
THE ALGORITHM 

Finally, we make a few comments about the accuracy 
and the performance of the algorithm. One indicator for 
the accuracy of the algorithm is the distance between 
the time-evolved state and the state with reduced vir- 
tual dimension. For the time-evolution of the Mott- 
distribution that was discussed in section fVTl this quan- 
tity is plotted in fig. [5J We find that the distance is typi- 
cally of order 10 -3 for D = 2 and of order 10~ 4 for D = 3 
and D = 4. Another quantity we monitor is the total 
number of particles (N) . Since this quantity is supposed 
to be conserved during the whole evolution, its fluctia- 



tions indicate the reliability of the algorithm. From the 
inset in fig. [9j the fluctuations of the particle number in 
case of the time-evolution of the Mott-distribution can 
be gathered. We find that these fluctuations are at most 
of order 10~ 5 . 

The main bottleneck for the performance of the algo- 
rithm is the scaling of the number of required multipli- 
cations with the virtual dimension D. As mentioned in 
section IIV1 the number of required multiplications is of 
order D 12 . Our simulations were run on a workstation 
with a 3.0 GHz Intel Xeon processor. On such a sys- 
tem, one evolution step on a 11 x 11-lattice with D = 5 
required a computing time of 55 hours. Another bottle- 
neck for the algorithm forms the scaling of the required 
memory with the virtual dimension D - which is of or- 
der D 8 . The simulation on a 11 x 11-lattice with D = 5 
thereby required a main memory of 2 GB. These bot- 
tlenecks make it difficult at the moment to go beyond 
a virtual dimension of D = 5. Nonetheless, a virtual 
dimension of D = 5 is expected to yield good results 
for many problems already. We intend to overcome the 
limitations of time and space by distributing tensor con- 
tractions among several processors in a future project. 

VIII. CONCLUSIONS 

Summing up, we have studied the system of hard- 
core bosons on a 2-D lattice using a variational method 
based on PEPS. We have thereby investigated the ground 
state properties of the system and its responses to sudden 
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changes in the parameters. We have compared our results 
to results based on the Gutzwiller ansatz. We have ob- 
served that the Gutzwiller ansatz predicts very well the 
density distribution of the particles. However, the mo- 
mentum distribution obtained from the Gutzwiller ansatz 
is, though qualitatively similar, quantitatively clearly dif- 
ferent from the distribution obtained from the PEPS 
ansatz. In addition, the PEPS and the Gutzwiller ansatz 



are very different in the prediction of time evolutions. 
We conclude that the Gutzwiller ansatz has to be ap- 
plied carefully in these cases. The simulations done in 
this paper give a clear demonstration of the power of 
the PEPS-approach, both for finding ground states in 
higher-dimensional quantum spin systems and for simu- 
lating real-time evolution. 
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